#
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
topocmap = plt.cm.gist_earth
topocmap.set_bad('w', alpha=0)
#First the greater area map
with nc(os.path.join(geodataf,'topo_2.nc')) as topof:
Gtopo=np.ma.masked_equal(topof.variables['topo'][:],0)
Gtopolon=topof.variables['longitude'][:]
Gtopolat=topof.variables['latitude'][:]
#Now only the tiwi-islands manp
with nc(os.path.join(geodataf,'Topo.nc')) as topof:
topo=np.ma.masked_equal(topof.variables['topo'][:],0)
topolon=topof.variables['longitude'][:]
topolat=topof.variables['latitude'][:]
with nc(CPOLF) as fnc:
lon=fnc.variables['longitude'][:]
lat=fnc.variables['latitude'][:]
fig = plt.figure()
ax = fig.add_subplot(111)
ls = LightSource(azdeg=315, altdeg=45)
M = Basemap(projection='rotpole', llcrnrlat=min(Gtopolat), llcrnrlon=min(Gtopolon), urcrnrlat=max(Gtopolat), urcrnrlon=max(Gtopolon),
resolution='c', area_thresh=1, ax=ax, lon_0=0, o_lon_p=0., o_lat_p=90.,)
borders = (([122.963, 139.055], [-17.091, -7.443]),
([127.385, 134.693], [-14.505, -10.005]),
([129.557, 132.529], [-13.757, -10.757]))
text=(('4km res.', 38), ('1.33 km res.', 32), ('444 m res.', 26))
M.pcolormesh(Gtopolon, Gtopolat, ls.hillshade(Gtopo[:], vert_exag=1), cmap= topocmap, ax=ax)
#M.drawmapscale(M.boundarylons[0]-0.7, M.boundarylats[0]+0.7, max(topolon), min(topolat), 10,
# barstyle='fancy', fontsize=19, labelstyle='simple')
for conf in range(3):
lons,lats = borders[conf]
txt, fts = text[conf]
x, y = M([lons[0], lons[0], lons[1],lons[1]], [lats[0], lats[1], lats[1], lats[0]])
xy = list(zip(x,y))
ttt = plt.text(*M(lons[0]+0.1, lats[0]+0.05), txt, fontsize=fts)
ttt.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='none'))
poly = Polygon(xy, alpha=1, edgecolor='k', lw=3, facecolor='none')
plt.gca().add_patch(poly)
#M.drawcoastlines()
M.scatter(*M(131.043101, -12.250323), marker='*',color='b', s=550)
plt.text(*M(borders[0][0][0]+0.15, borders[0][-1][1]-0.4), 'b)', fontsize=32)
axins = zoomed_inset_axes(ax, 7, loc=2)
axins.set_xlim(min(topolon), max(topolon))
axins.set_ylim(min(topolat), max(topolat))
tmap = Basemap(projection='cyl',llcrnrlat=min(lat), llcrnrlon=min(lon), urcrnrlat=max(lat), urcrnrlon=max(lon),
resolution='c', area_thresh=1, ax=axins)
tmap.pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap= topocmap)
mark_inset(ax, axins, loc1=3, loc2=4, fc="none", ec='k')
plt.annotate('a)', xy=(0.01 ,0.95), xycoords='axes fraction', fontsize=32)
lons, lats = [min(Gtopolon), max(Gtopolon)], [min(Gtopolat), max(Gtopolat)]
axins2 = inset_axes(ax, width='15%', height='15%', loc=3)
Worldmap= Basemap(projection='ortho', lat_0=Gtopolat[int(len(Gtopolat)/2)], lon_0=Gtopolon[int(len(Gtopolon)/2)],
ax=axins2, height=240*10**3)
Worldmap.drawcoastlines(color='white')
Worldmap.fillcontinents(color='gray')
bx, by = Worldmap(M.boundarylons, M.boundarylats)
xy = list(zip(bx,by))
poly = Polygon(xy, alpha=1, edgecolor='none', lw=0, facecolor='darkcyan')
Worldmap.ax.add_patch(poly)
warning: width and height keywords ignored for Orthographic projection
<matplotlib.patches.Polygon at 0x7f4c55b05278>
To investigate the links between Monsoon stage and Rainfall extremes in the considered area it is worthwhile having a look at the general climatology of rainfall and investigate the importance of the diurnal cycle in the stduy domain.
The relative importance of the diurnal cycle is investigated in two ways:
Because the cpol data contains holes (missing data) we chose to create a Lomb-Scargle periodigram rather than a a fourier decomposition of the data. Lomb-Scargle periodigrams are created by least-spaure fitting sinosoids to the data. The LombScargle Class provided by astropy's stats module will create the periodogram. For best practice we create a function that calculates the fraction of the spectral variance.
#Plot the output map
fig = plt.figure()
#ax1 = fig.add_subplot(111)
ax1 = fig.add_subplot(131)
ax1.set_title('rainfall climatology (1998 - 2017)', size=28)
im = m.pcolormesh(lon, lat, mask*np.nanmean(cpol24h,axis=0)*24, vmin=0, cmap=colmap.GMT_drywet, ax=ax1)
m.drawcoastlines()
cbar = m.colorbar(im,location='bottom',pad="5%")
cbar.set_label('Avg. rain-rate [mm/d]', size=28)
cbar.ax.tick_params(labelsize=24)
#'''
ax2 = fig.add_subplot(132)
ax2.semilogx(f, psd * 9, lw=4)
ax2.set_xlim(0.25,60)
ax2.set_ylim(0,0.6)
ax2.set_xlabel('Periods ($\\tau$) [d]',size=28)
ax2.set_ylabel('Spectral variance [%]',size=28)
ax2.tick_params(labelsize=28)
ax3 = fig.add_subplot(133)
im2 = m.pcolormesh(lon, lat, np.ma.masked_invalid(spec_map)*2 , vmin=0, vmax=0.8,cmap='Blues', ax=ax3)
m.drawcoastlines()
cbar = m.colorbar(im2,location='bottom',pad="5%")
cbar.set_label('Spectral variance of $\\tau$ = 1 d [%]', size=28)
cbar.ax.tick_params(labelsize=24)
# %%
fig.subplots_adjust(right=0.99, bottom=0.83, top=0.99,left=0.01, hspace=0.05, wspace=0.1)
#'''
#Plot the diurnal Cycle
fig, ax = plt.subplots(1, figsize=(12,6))
ax.set_ylabel('Rain-rate mm/d', fontsize=22)
ax.set_xlabel('Local-time', fontsize=22)
ax.plot(idx, dc_burst.values, 'fuchsia', linewidth=2.0, label='Bursts')
ax.plot(idx, dc_break.values, 'g', linewidth=2.0, label='Breaks')
ax.plot(idx, (dc_break+dc_burst)/2, 'firebrick', linewidth=2.0, label='All')
ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%_H'))
ax.legend(loc=0)
ax.grid(linewidth=0)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
The above plots show that diurnal cycle is the most important mode of variability in the considered area. The island thunderstrom hector is the most extreme convective system that occurs at the time-scale at perdiods of 24 h. The Hector 'signal' is most pronounced during the Monsoon break period. Before studying Hector we split the CPOL area into break and burst period and investigate the occurrence of extremes during Monsoon break and burst.
First we define extreme events by applying the 99th percentile threshold of the total 10min dataset. Then we define the break and burst events by the definition of Narsey et al. and calculate the number of burst and the number of break events during the CPOL era
# Now Plot the maps
maxval =[(0,15), (0,15), (-12,12), (0, 150), (0, 150), (-200, 200)]
names=['Avg Rain during breaks', 'Avg Rain during bursts', 'Difference (breaks - bursts)',
'Extremes during breaks', 'Extremes during bursts', 'Difference (breaks - bursts)']
cmap_list = [(colmap.GMT_drywet,'k'), (colmap.GMT_drywet,'k'), (colmap.GMT_polar_r,'k'),
(mpl.cm.magma,'w'), (mpl.cm.magma,'w'), (colmap.GMT_polar_r,'k')]
C = [Mean[1]*1.25, Mean[2], Mean[1]*1.25 - Mean[2], D[1], D[2], D[1] - D[2]]
fig = plt.figure()
mask = ((D[0] * 0 )+1)
im=[]
for i, idx in enumerate(C):
ax = fig.add_subplot(2,3,i+1)
data = np.zeros([len(lat),len(lon)])
colm, coast = cmap_list[i]
try:
im.append(m.pcolormesh(lon, lat, mask*C[i].filled(0),vmin=maxval[i][0],vmax=maxval[i][1],cmap=colm))
except AttributeError:
im.append(m.pcolormesh(lon, lat, mask*C[i],vmin=maxval[i][0],vmax=maxval[i][1], cmap=colm))
m.drawcoastlines(color=coast)
#cbar=m.colorbar(im,location='bottom',pad='5%')
#cbar.set_label(u'Freq. of events above 99th% [%]',size=18)
#cbar.ax.tick_params(labelsize=18)
ax.set_title('%s'%names[i],size=28)
fig.subplots_adjust(right=0.95, bottom=.45, top=0.95,left=0.05, hspace=0, wspace=0)
cbar_ax1 = fig.add_axes([0.05, 0.72,0.58, 0.01])
cbar1 = fig.colorbar(im[0], cax=cbar_ax1, orientation='horizontal')
cbar1.ax.tick_params(labelsize=24)
cbar1.set_label(u'Rain-rate [mm/d]',size=28)
cbar_ax2 = fig.add_axes([0.67, 0.72,0.27, 0.01])
cbar2 = fig.colorbar(im[2], cax=cbar_ax2, orientation='horizontal')
cbar2.ax.tick_params(labelsize=24)
cbar2.set_label(u'Difference',size=28)
cbar_ax3 = fig.add_axes([0.05, 0.47,0.58, 0.01])
cbar3 = fig.colorbar(im[3], cax=cbar_ax3, orientation='horizontal')
cbar3.ax.tick_params(labelsize=24)
cbar3.set_label(u'N. of events above 99th%',size=28)
cbar_ax4 = fig.add_axes([0.67, 0.47,0.27, 0.01])
cbar4 = fig.colorbar(im[-1], cax=cbar_ax4, orientation='horizontal')
cbar4.ax.tick_params(labelsize=24)
cbar4.set_label(u'Difference',size=28)
To investigate how well Extreme Events caused by local tropicsl Island Thunderstorms are represented in a numerical model we choose a period of 8 consecutive occurring Hector events over the Tiwi Island during the Monsoon break period in November 2006. To increase the sample size for more robust statistics we create ensemble of 8 members. The ensemble is created by 8 different initialisation times.
Now we compare the ensemble timeseries of the simulated rainfall for the 1.33km and 0.44km Simulations with the CPOL observations.
#Plot area avg TS
import matplotlib.dates as mdates
myFmt = mdates.DateFormatter('%d/%m %H LT')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(20,10), dpi=72)
ax = fig.add_subplot(111)
obs_time_1 = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
obs_time_2 = pd.DatetimeIndex(OBS.dataset[2].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
um044_time = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
um133_time = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
## First plot the area avg of the cpol data
#ax.plot(OBS.dataset[0].coords['t'], np.nanmean(OBS.dataset[0].variables['lsrain'].values[:,0,:]*mask.filled(np.nan),axis=(1,2)),
# label='CPOL', color='firebrick')
#ax.plot(obs_time_2, OBS.dataset[2].variables['lsrain'].mean(axis=(1,2)), label='Cmorph 8km',
# color='fuchsia', linestyle='-')
ax.plot(obs_time_1, OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)), linestyle='-',
color='lightseagreen', label='CPOL')
#############
ax.plot(um133_time, UM133ens.variables['lsrain'].mean(axis=(0,2,3,4)), color='fuchsia',
label='UM 1.33km ensmean', linestyle='--')
#ax.fill_between(UM133ens.coords['t'].values,
# UM133ens.variables['lsrain'].mean(axis=(2,3,4)).min(axis=0).values,
# UM133ens.variables['lsrain'].mean(axis=(2,3,4)).max(axis=0).values, color='fuchsia',
# alpha=0.2)
############
ax.plot(um044_time, UM044ens.variables['lsrain'].mean(axis=(0,2,3,4)), color='fuchsia',
label='UM 0.44km ensmean')
ax.fill_between(um044_time,
UM044ens.variables['lsrain'].mean(axis=(2,3,4)).min(axis=0).values,
UM044ens.variables['lsrain'].mean(axis=(2,3,4)).max(axis=0).values, color='cornflowerblue',
alpha=0.2)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
ax.set_ylim(0.0,0.42)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=18)
ax.xaxis.set_major_formatter(myFmt)
ax.grid(color='r', linestyle='-', linewidth=0)
Now lets look at some animations of UM rainfall occuring over the Tiwi-islands and compare it to the cpol observations
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Ens-2.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
#Create the avg. maps of rainfall 0.44km
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.35, top=0.95,left=0.05, hspace=0, wspace=0)
UM = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m[0].drawcoastlines()
cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['Obs']
ax[-1].set_title('Obs', fontsize=24)
for i in range(1,9):
ax.append(fig.add_subplot(3,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=100, ax=ax[-1]))
m[-1].drawcoastlines()
#m[-1].shadedrelief()
im.append(m[-1].pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append(datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime('%e %b %H UTC'))
ax[-1].set_title(tit[-1],fontsize=24)
#Create the avg maps of rainfall for 1.33km
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.35, top=0.95,left=0.05, hspace=0, wspace=0)
UM = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m[0].drawcoastlines()
cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['Obs']
ax[-1].set_title('Obs', fontsize=24)
for i in range(1,9):
ax.append(fig.add_subplot(3,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=100, ax=ax[-1]))
m[-1].drawcoastlines()
#m[-1].shadedrelief()
im.append(m[-1].pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append(datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime('%e %b %H UTC'))
ax[-1].set_title(tit[-1],fontsize=24)
#Plet the avg maps for comparison
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.01, top=0.6,left=0.01, hspace=0, wspace=0)
UM_2 = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
UM_1 = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01)
data = [obs_data[1:-1]/6,
np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
None,
np.ma.masked_less(UM_1.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0],
np.ma.masked_less(UM_2.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0]
]
data.append(None)
data.append(data[0] - data[1])
data.append(data[0] - data[2])
UM_1min = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_2min = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_1max = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
UM_2max = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
#data.append(None)
#data.append(UM_1max - UM_1min)
#data.append(UM_1max - UM_1min)
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM_1.coords['lat'][:], UM_1.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m[0].drawcoastlines()
#cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
#cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
#cbar.ax.tick_params(labelsize=24)
#cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['CPOL avg.', 'UM 0.44km ens avg.', 'UM 1.33km ens avg.',
None, 'UM 0.44 ens std.', 'UM 1.33km ens std.',
None, 'CPOL - UM 0.44 ens avg.' ,'CPOL - UM 1.33km ens avg.',
None, 'UM 0.44 ens range', 'UM 1.33km ens range']
colors = {1:colm, 2:colm, 4:colm, 5:colm,
7:colmap.GMT_polar_r, 8:colmap.GMT_polar_r,
10:colm, 11:colm}
Range = {1:(0.0, 25/6), 2:(0.0, 25/6), 4:(0,8/6), 5:(0,8/6), 7:(-20/6,20/6), 8:(-20/6,20/6), 10:(0,30/6), 11:(0,30/6)}
height = {2:0.43, 5:0.23, 8:0.03, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=24)
cbar_ax = [fig.add_axes([0.95, height[2], 0.01, 0.15])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=24)
cbar[-1].set_label('Rain-rate [mm/d]',size=24)
for i in range(1,len(data)):
if data[i] is not None:
ax.append(fig.add_subplot(3,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=100, ax=ax[-1]))
m[-1].drawcoastlines()
#m[-1].shadedrelief()
im.append(m[-1].pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
ax[-1].set_title(tit[i],fontsize=24)
if i in (5, 8, 11):
cbar_ax.append(fig.add_axes([0.95, height[i], 0.01, 0.15]))
cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
cbar[-1].ax.tick_params(labelsize=24)
cbar[-1].set_label('Rain-rate [mm/d]',size=24)
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(20,10), dpi=72)
ax = fig.add_subplot(111)
####################
cpol_1 = OBS.dataset[1].groupby('t.hour').mean(dim='t').mean(dim=('y','x'))
obs_S = pd.Series(cpol_1.variables['lsrain'].values, index=(cpol_1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
ax.plot(obs_S.index, obs_S.values, linestyle='-', color='lightseagreen', label='CPOL')
##################
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).mean(dim='ens')
um_2 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).min(dim='ens')
um_3 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).max(dim='ens')
um_044 = pd.DataFrame({'mean':um_1.variables['lsrain'].values,
'min' :um_2.variables['lsrain'].values,
'max' :um_3.variables['lsrain'].values},
index=(um_1.variables['hour'].values + 9) % 24).sort_index(inplace=False)
ax.plot(um_044.index, um_044['mean'].values, color='fuchsia', linestyle='-', label='UM 0.44km ensmean')
###################
um_4 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).mean(dim='ens')
um_5 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).min(dim='ens')
um_6 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).max(dim='ens')
um_133 = pd.DataFrame({'mean':um_4.variables['lsrain'].values,
'min' :um_5.variables['lsrain'].values,
'max' :um_6.variables['lsrain'].values},
index=(um_4.variables['hour'].values + 9) % 24).sort_index(inplace=False)
ax.plot(um_133['mean'].index, um_133['mean'].values, color='fuchsia', linestyle='--', label='UM 1.33km ensmean')
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
#####################
ax.fill_between(um_044.index, um_044['min'], um_044['max'], color='cornflowerblue', alpha=0.2, label='UM 0.44km ens. spread')
ax.fill_between(um_133.index, um_133['min'], um_133['max'], color='coral', alpha=0.2, label='UM 1.33km ens. spread')
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=18)
ax.set_xlabel('Local time [h]', fontsize=18)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
ax.grid(color='r', linestyle='-', linewidth=0)
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Diurnal-2.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
#Plot the maps
plt.close()
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
um_2 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
obs = OBS.dataset[1].groupby('t.hour').mean(dim='t')
time = (um_1.hour + 9) % 24
data = [obs, um_2, um_1]
#fig = plt.figure(figsize=(25,15))
lat1, lon1 = OBS.dataset[1].variables['latitude'].values[:,0], OBS.dataset[1].variables['longitude'].values[0,:]
lat2, lon2 = um_1.coords['lat'][:], um_1.coords['lon'][:]
#fig.subplots_adjust(right=0.94, bottom=0.025, top=0.67,left=0.01, hspace=0, wspace=0)
first = True
val_max = 3
colm = colmap2.Blues
colm.set_under('w', alpha=0)
m, ax, im = [], [], []
for i in range(24):
break
fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
if first:
for col in range(3*2):
if col != 3:
ax.append(fig.add_subplot(2,3,col+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1, ax=ax[-1]))
m[-1].drawcoastlines()
###MEAN
ax[0].set_title('CPOL',fontsize=22)
ax[1].set_title('UM 1.33km mean', fontsize=22)
ax[2].set_title('UM 0.44km mean', fontsize=22)
ax[-2].set_title('UM 1.33km std', fontsize=22)
ax[-1].set_title('UM 0.44km std', fontsize=22)
im.append(m[0].pcolormesh(lon1, lat1, obs['lsrain'].values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
im.append(m[1].pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
im.append(m[2].pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
###STD
im.append(m[-2].pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
im.append(m[-1].pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
first = False
cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
else:
###MEAN
im[0].set_array(obs['lsrain'].values[i].ravel())
im[1].set_array(um_2['lsrain'].mean(dim='ens').values[i].ravel())
im[2].set_array(um_1['lsrain'].mean(dim='ens').values[i].ravel())
###STD
im[-2].set_array(um_2['lsrain'].std(dim='ens').values[i].ravel())
im[-1].set_array(um_1['lsrain'].std(dim='ens').values[i].ravel())
cbar.set_label('Rain-rate at %2i LT [mm/h]'%((i+9)%24),size=24)
fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
#break
day = (time[1:5] - 9) % 24
night = (time[5:5+4] - 9) % 24
obs_diff = obs['lsrain'][day].sum(dim='hour') - obs['lsrain'][night].sum('hour')
um_1_diff = um_1['lsrain'][day].sum(dim='hour') - um_1['lsrain'][night].sum('hour')
um_2_diff = um_2['lsrain'][day].sum(dim='hour') - um_2['lsrain'][night].sum('hour')
mask = np.ma.masked_invalid(obs['lsrain'].values) * 0 + 1
obs_max = (((np.argmax(obs['lsrain'].values * mask, axis=0)*mask[0] + 9 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9) % 24 / 2).round(0) * 2
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9) % 24 / 2).round(0) * 2
colm = CubeYF_20.get_mpl_colormap(N=24, gamma=2.0)
#colm = qualitative.Paired[12].get_mpl_colormap(N=24, gamma=2.0)
m, ax = [], []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km ens. mean', 'UM 0.44km ens. mean']
fig = plt.figure(figsize=(25,15))
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.67,left=0.01, hspace=0, wspace=0)
cbar_ax = fig.add_axes([0.12, 0.17, 0.74, 0.02])
for i in range(3):
ax.append(fig.add_subplot(1,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1, ax=ax[-1]))
m[-1].drawcoastlines()
try:
im = m[-1].pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
except TypeError:
im = m[-1].pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
ax[-1].set_title(names[i], fontsize=22)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Local Time [h]', size=24)
#Plot the pdfs
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.UM044['shape']*3.4], loc = np.log(LogNorm.UM044['scale']*1.9))#, shape # mu, sigma
med = np.median(um044vec.loc[um044vec>0].values)
X = np.linspace(0,65,200)
#b = 2.5
nbins = 100
def fitpdf(x, b):
a=med*100
return ((b/a) * (x/a)**(b-1)) / (1+(x/a)**b)**2
#popt, pcov = curve_fit(powerlaw, um044cnt['Rain-rate'].values[1:],
# um044cnt['count'].values[1:])
popt, pcov = PowerFit.UM044['popt']/4, PowerFit.UM044['pcov']
fig = plt.figure()
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um044vec.loc[um044vec>0], um133vec.loc[um133vec>0], obs_vec.loc[obs_vec>0])
facecolors = ['fuchsia', 'firebrick', 'lightseagreen']
#n1_1, bins1_1, patches1_1 = ax2.hist(, bins=nbins, normed=1, facecolor='firebrick', alpha=0.25)
#n2, bins2, patches2 = ax2.hist(, bins=nbins, normed=1, facecolor='lightseagreen', alpha=0.25)
ax2.set_ylim(0,1.4)
ax2.set_yticks([])
ax.plot(np.linspace(0.01, 100,200), fitpdf(np.linspace(0.01,100,200),popt[0]), lw=4, color='navy',
label='Log-logistic $\\beta = %02.2e$'%popt[0])
ax.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
ax.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
ax.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#ax.scatter(um044cnt['Rain-rate'].values, um044cnt['count'].values, c='red', s=125, alpha=0.24)
#ax.scatter(obscnt['Rain-rate'].values, obscnt['count'].values, c='purple', s=125, alpha=0.24)
#ax.plot(np.linspace(0, 50,200), um044pdf, lw=4, label='UM 0.44km ens', color='lightseagreen')
#ax.plot(np.linspace(0, 50,200), obspdf, lw=4, label='CPOL', color='fuchsia')
ax.plot(np.linspace(0.0,50,200), dist.pdf(np.linspace(0.0, 50,200)), color='firebrick', alpha=0.8,
lw=4, label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
n1, bins1, patches1 = ax2.hist(histdata, bins=nbins, color=facecolors, normed=1, alpha=0.25)
#fitted_pl.plot_pdf(ax=ax,color='k',lw=4)
#ax.plot(np.linspace(0.1, 100), um044dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 0.44km area-avg')
#ax.plot(np.linspace(0.1, 100), um133dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 1.33km area-avg')
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.legend(loc=4, fontsize=28)
ax.tick_params(labelsize=28)
ax.set_ylim(0.,.4)
ax.set_xlim(0.,40)
ax.set_ylabel('Density', fontsize=28)
ax.set_xlabel('Rain-rate mm/h', fontsize=28)
axin = inset_axes(ax, width = "60%", height = "60%", loc=1)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#axin.plot(X, um044pdf,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#axin.plot(X, obspdf,'-', label='CPOL',lw=4, color='lightseagreen')
axin.plot(np.linspace(0.01,80,200), dist.pdf(np.linspace(0.01, 80,200)), lw=4, color='firebrick', alpha=0.8,
label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=24)
axin.set_xlim(0.01, 90)
axin.set_ylim(0.000025, 1)
axin.grid(color='r', linestyle='-', linewidth=0)
axin.annotate("b)", (0.01,0.03), xycoords="axes fraction", fontsize=28)
ax.grid(color='r', linestyle='-', linewidth=0)
axin.plot(np.linspace(0.01, 50,200), fitpdf(np.linspace(0.01,50,200),popt[0]), lw=4, color='navy')
axin.set_ylabel('$\\log$(Density)', fontsize=24)
axin.set_xlabel('$\\log$(Rain-rate)', fontsize=24)
<matplotlib.text.Text at 0x7f4c403015f8>
Now get the percentiles for the ensemble simulation and compare it against the observations
#Plot the percentages
fig=plt.figure()
ax = fig.add_subplot(111)
ax.plot(PERC.index, PERC['Obs'].values, color='lightseagreen', linestyle='-', label='CPOL',lw=3)
ax.plot(PERC.index, PERC['UM 1.33km'].values, color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], color='lightseagreen', linestyle='-', label='CPOL',lw=3)
axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=24)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])
axin.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
#axin.get_xaxis().get_major_formatter().labelOnlyBase = False
axin.grid(color='r', linestyle='-', linewidth=0)
axin.annotate("b)", (0.01,0.03), xycoords="axes fraction", fontsize=28)
axin.set_xticks(range(50,100,10), minor=False)
ax.grid(color='r', linestyle='-', linewidth=0)
axin.set_ylabel('Rain-rate[mm/h]', fontsize=24)
axin.set_xlabel('Percentile [%]', fontsize=24)
<matplotlib.text.Text at 0x7f4c38acc7b8>
##Intermezzo
first = True
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.05, top=0.95,left=0.05, hspace=0, wspace=0)
data1 = np.ma.masked_less(OBS.dataset[-1].variables['lsrain'][:].values, 0.02)
data2 = np.ma.masked_less(OBS.dataset[1].variables['lsrain'][:].values, 0.02)
t1 = OBS.dataset[-1].variables['t'][:]
t2 = OBS.dataset[1].variables['t'][:]
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
#for tt, tstep in enumerate(tsteps):
tt = 0
fname = os.path.join(outdir, 'Comparison_%05i.png'%tt)
#sys.stdout.flush()
#sys.stdout.write('\r Creating %s .... '%os.path.basename(fname))
#sys.stdout.flush()
lon1, lat1 = OBS.dataset[-1].variables['longitude'][0,:], OBS.dataset[-1].variables['latitude'][:,0]
lon2, lat2 = OBS.dataset[1].variables['longitude'][0,:], OBS.dataset[1].variables['latitude'][:,0]
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1), resolution='i',
area_thresh=1, ax=ax1)
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='i',
area_thresh=1, ax=ax2)
im1 = m1.pcolormesh(lon1, lat1, data1[tt].filled(-1),vmin=0.0,vmax=30,cmap='Blues')
im2 = m2.pcolormesh(lon2, lat2, data2[tt].filled(-1),vmin=0.0,vmax=30,cmap='Blues')
m1.drawcoastlines()
m2.drawcoastlines()
cbar_ax = fig.add_axes([0.14, 0.35, 0.74, 0.01])
cbar=fig.colorbar(im1, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Rain-rate [mm/h]',size=24)
ax2.set_title('Old CPOL version at %s LT' %tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
ax1.set_title('New CPOL version at %s LT'%tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
def ani(tt):
im1.set_array((data1[tt]).filled(-1).ravel())
im2.set_array((data2[tt]).filled(-1).ravel())
ax2.set_title('Old CPOL version at %s LT' %tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
ax1.set_title('New CPOL version at %s LT'%tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
return(im1, im2)
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, ani, frames=100, interval=20, blit=True)
#fig.savefig(fname, bbox_inches='tight', format='png', dpi=72)
#break
#sys.stdout.write('\n')
#Visualize intermezzo
import io
import base64
from IPython.display import HTML
viddir = os.path.join(os.getenv('HOME'),'Data/Extremes/CPOL/Plot/Versions')
outvid=os.path.join(viddir,'WeekOfHector.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
Now we compare the distributions of percentiles during break and burst periods. We also add
#Plot percentiles for breaks, bursts and all
colors = {'10min':'firebrick'}# '1h':'lightseagreen', '3h':'fuchsia' , '6h':'cornflowerblue'}
fig=plt.figure()
ax = fig.add_subplot(111)
for con in (('break','-'), ('burst','--'), ('all','-.')):
for gr in colors.keys():
df = pd.read_hdf(PercF, con[0])[gr]
ax.plot(df.index, df.values, color=colors[gr], linestyle=con[-1], label='%s (%s)'%(gr, con[0]))
ax.set_xlim(1,99.9999)
#ax.set_ylim(0.11,)
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_xlabel('Percentile [%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)